rm(list=ls())

set.seed(1234)

library(plyr)
library(ggplot2)
library(car)
library(psych)
library(texreg)
library(effectsize)
library(xtable)

setwd(".../Data/Analysis")

#Read in individual level datasets 

Mex<-read.csv("MexAgg.csv")
Peru<-read.csv("PeruAgg.csv")
UK<-read.csv("UKAgg.csv")
Spain<-read.csv("SpainAgg.csv")
Port<-read.csv("PortAgg.csv")
USA<-read.csv("USAgg.csv")
Brazil<-read.csv("BrazilAgg.csv")
Arg<-read.csv("ArgAgg.csv")
NZ<-read.csv("NZAgg.csv")
NorAusFrance<-read.csv("NorAusFranceAgg.csv")
Nor<-subset(NorAusFrance, QCOUNTRY == 'Norway')
France <-subset(NorAusFrance, QCOUNTRY == 'France')
Aus <-subset(NorAusFrance, QCOUNTRY == 'Australia')

#Assess quota penality for each country 

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.1) # move them .05 to the left and right

breaks<-seq(1,5,1)


## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

#Mexico
Mex <-subset(Mex, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Mex <- summarySE(Mex, measurevar="SubLegStand", groupvars=c("treat"))
MexDiffSub<-data_sum_Mex$SubLegStand[3] - data_sum_Mex$SubLegStand[2]
MexSESub<-data_sum_Mex$se[3]
Mexico<-c("Mexico", MexDiffSub, MexSESub, 50) #ADD quota threshold from Table 1

#Peru
Peru <-subset(Peru, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Peru <- summarySE(Peru, measurevar="SubLegStand", groupvars=c("treat"))
PeruDiffSub<-data_sum_Peru $SubLegStand[3] - data_sum_Peru$SubLegStand[2]
PeruSESub<-data_sum_Peru$se[3]
Peru <-c("Peru", PeruDiffSub, PeruSESub, 50) #ADD quota threshold from Table 1

#UK
UK <-subset(UK, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_UK <- summarySE(UK, measurevar="SubLegStand", groupvars=c("Vignettes_DO"))
UKDiffSub<-data_sum_UK$SubLegStand[3] - data_sum_UK$SubLegStand[2]
UKSESub<-data_sum_UK$se[3]
UK <-c("UK", UKDiffSub, UKSESub, 16) #ADD quota threshold from Table 1

#Spain
Spain <-subset(Spain, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Spain <- summarySE(Spain, measurevar="SubLegStand", groupvars=c("Vignettes_DO"))
SpainDiffSub<-data_sum_Spain $SubLegStand[3] - data_sum_Spain$SubLegStand[2]
SpainSESub<-data_sum_Spain$se[3]
Spain <-c("Spain", SpainDiffSub, SpainSESub, 40) #ADD quota threshold from Table 1

#USA
USA <-subset(USA, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_USA <- summarySE(USA, measurevar="SubLegStand", groupvars=c("treat"))
USADiffSub<-data_sum_USA$SubLegStand[3] - data_sum_USA$SubLegStand[2]
USASESub<-data_sum_USA$se[3]
USA <-c("USA", USADiffSub, USASESub, 0) #ADD quota threshold from Table 1

#Portugal
Port <-subset(Port, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Port <- summarySE(Port, measurevar="SubLegStand", groupvars=c("Vignettes_DO"))
PortDiffSub<-data_sum_Port$SubLegStand[3] - data_sum_Port$SubLegStand[2]
PortSESub<-data_sum_Port$se[3]
Port <-c("Portugal", PortDiffSub, PortSESub, 33) #ADD quota threshold from Table 1

#Norway
Norway <-subset(Nor, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Norway <- summarySE(Norway, measurevar="SubLegStand", groupvars=c("treat"))
NorDiffSub<-data_sum_Norway$SubLegStand[3] - data_sum_Norway$SubLegStand[2]
NorSESub<-data_sum_Norway$se[3]
Norway <-c("Norway", NorDiffSub, NorSESub, 32) #ADD quota threshold from Table 1

#Australia
Aus <-subset(Aus, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Aus <- summarySE(Aus, measurevar="SubLegStand", groupvars=c("treat"))
AusDiffSub<-data_sum_Aus$SubLegStand[3] - data_sum_Aus$SubLegStand[2]
AusSESub<-data_sum_Aus$se[3]
Aus <-c("Australia", AusDiffSub, AusSESub, 23) #ADD quota threshold from Table 1

#France
France <-subset(France, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_France <- summarySE(France, measurevar="SubLegStand", groupvars=c("treat"))
FranceDiffSub<-data_sum_France$SubLegStand[3] - data_sum_France$SubLegStand[2]
FranceSESub<-data_sum_France$se[3]
France <-c("France", FranceDiffSub, FranceSESub, 50) #ADD quota threshold from Table 1

#Argentina
Argentina <-subset(Arg, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Arg <- summarySE(Argentina, measurevar="SubLegStand", groupvars=c("treat"))
ArgDiffSub<-data_sum_Arg$SubLegStand[3] - data_sum_Arg$SubLegStand[2]
ArgSESub<-data_sum_Arg$se[3]
Argentina <-c("Argentina", ArgDiffSub, ArgSESub, 50) #ADD quota threshold from Table 1

#Brazil
Brazil <-subset(Brazil, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Brazil <- summarySE(Brazil, measurevar="SubLegStand", groupvars=c("treat"))
BrazilDiffSub<-data_sum_Brazil$SubLegStand[3] - data_sum_Brazil$SubLegStand[2]
BrazilSESub<-data_sum_Brazil$se[3]
Brazil <-c("Brazil", BrazilDiffSub, BrazilSESub, 30) #ADD quota threshold from Table 1

#NZ
NZ <-subset(NZ, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_NZ <- summarySE(NZ, measurevar="SubLegStand", groupvars=c("treat"))
NZDiffSub<-data_sum_NZ$SubLegStand[3] - data_sum_NZ$SubLegStand[2]
NZSESub<-data_sum_NZ$se[3]
NZ <-c("New Zealand", NZDiffSub, NZSESub, 25) #ADD quota threshold from Table 1

penaltySub<-as.data.frame(rbind(Mexico, Peru, UK, Spain, USA, Port, Norway, Aus, France, Argentina, Brazil, NZ))

##Do the same procedure for procedural legitimacy

#re-read in data

Mex<-read.csv("MexAgg.csv")
Peru<-read.csv("PeruAgg.csv")
UK<-read.csv("UKAgg.csv")
Spain<-read.csv("SpainAgg.csv")
Port<-read.csv("PortAgg.csv")
USA<-read.csv("USAgg.csv")
Brazil<-read.csv("BrazilAgg.csv")
Arg<-read.csv("ArgAgg.csv")
NZ<-read.csv("NZAgg.csv")
NorAusFrance<-read.csv("NorAusFranceAgg.csv")
Nor<-subset(NorAusFrance, QCOUNTRY == 'Norway')
France <-subset(NorAusFrance, QCOUNTRY == 'France')
Aus <-subset(NorAusFrance, QCOUNTRY == 'Australia')

#Mexico
Mex <-subset(Mex, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Mex <- summarySE(Mex, measurevar="ProLegStand", groupvars=c("treat"))
MexDiffSub<-data_sum_Mex$ProLegStand[3] - data_sum_Mex$ProLegStand[2]
MexSESub<-data_sum_Mex$se[3]
Mexico<-c("Mexico", MexDiffSub, MexSESub, 50) #ADD quota threshold from Table 1

#Peru
Peru <-subset(Peru, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Peru <- summarySE(Peru, measurevar="ProLegStand", groupvars=c("treat"))
PeruDiffSub<-data_sum_Peru $ProLegStand[3] - data_sum_Peru$ProLegStand[2]
PeruSESub<-data_sum_Peru$se[3]
Peru <-c("Peru", PeruDiffSub, PeruSESub, 50) #ADD quota threshold from Table 1

#UK
UK <-subset(UK, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_UK <- summarySE(UK, measurevar="ProLegStand", groupvars=c("Vignettes_DO"))
UKDiffSub<-data_sum_UK$ProLegStand[3] - data_sum_UK$ProLegStand[2]
UKSESub<-data_sum_UK$se[3]
UK <-c("UK", UKDiffSub, UKSESub, 16) #ADD quota threshold from Table 1

#Spain
Spain <-subset(Spain, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Spain <- summarySE(Spain, measurevar="ProLegStand", groupvars=c("Vignettes_DO"))
SpainDiffSub<-data_sum_Spain $ProLegStand[3] - data_sum_Spain$ProLegStand[2]
SpainSESub<-data_sum_Spain$se[3]
Spain <-c("Spain", SpainDiffSub, SpainSESub, 40) #ADD quota threshold from Table 1

#USA
USA <-subset(USA, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_USA <- summarySE(USA, measurevar="ProLegStand", groupvars=c("treat"))
USADiffSub<-data_sum_USA$ProLegStand[3] - data_sum_USA$ProLegStand[2]
USASESub<-data_sum_USA$se[3]
USA <-c("USA", USADiffSub, USASESub, 0) #ADD quota threshold from Table 1

#Portugal
Port <-subset(Port, Vignettes_DO == 'treat1' | Vignettes_DO == 'treat2' | Vignettes_DO == 'treat3') #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Port <- summarySE(Port, measurevar="ProLegStand", groupvars=c("Vignettes_DO"))
PortDiffSub<-data_sum_Port$ProLegStand[3] - data_sum_Port$ProLegStand[2]
PortSESub<-data_sum_Port$se[3]
Port <-c("Portugal", PortDiffSub, PortSESub, 33) #ADD quota threshold from Table 1

#Norway
Norway <-subset(Nor, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Norway <- summarySE(Norway, measurevar="ProLegStand", groupvars=c("treat"))
NorDiffSub<-data_sum_Norway$ProLegStand[3] - data_sum_Norway$ProLegStand[2]
NorSESub<-data_sum_Norway$se[3]
Norway <-c("Norway", NorDiffSub, NorSESub, 32) #ADD quota threshold from Table 1

#Australia
Aus <-subset(Aus, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Aus <- summarySE(Aus, measurevar="ProLegStand", groupvars=c("treat"))
AusDiffSub<-data_sum_Aus$ProLegStand[3] - data_sum_Aus$ProLegStand[2]
AusSESub<-data_sum_Aus$se[3]
Aus <-c("Australia", AusDiffSub, AusSESub, 23) #ADD quota threshold from Table 1

#France
France <-subset(France, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_France <- summarySE(France, measurevar="ProLegStand", groupvars=c("treat"))
FranceDiffSub<-data_sum_France$ProLegStand[3] - data_sum_France$ProLegStand[2]
FranceSESub<-data_sum_France$se[3]
France <-c("France", FranceDiffSub, FranceSESub, 50) #ADD quota threshold from Table 1

#Argentina
Argentina <-subset(Arg, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Arg <- summarySE(Argentina, measurevar="ProLegStand", groupvars=c("treat"))
ArgDiffSub<-data_sum_Arg$ProLegStand[3] - data_sum_Arg$ProLegStand[2]
ArgSESub<-data_sum_Arg$se[3]
Argentina <-c("Argentina", ArgDiffSub, ArgSESub, 50) #ADD quota threshold from Table 1

#Brazil
Brazil <-subset(Brazil, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_Brazil <- summarySE(Brazil, measurevar="ProLegStand", groupvars=c("treat"))
BrazilDiffSub<-data_sum_Brazil$ProLegStand[3] - data_sum_Brazil$ProLegStand[2]
BrazilSESub<-data_sum_Brazil$se[3]
Brazil <-c("Brazil", BrazilDiffSub, BrazilSESub, 30) #ADD quota threshold from Table 1

#NZ
NZ <-subset(NZ, treat == 1 | treat == 2 | treat == 3) #quota penalty is treatment 3 (quota) - treatment 2 (gender-balanced); negative value is the "penalty"
data_sum_NZ <- summarySE(NZ, measurevar="ProLegStand", groupvars=c("treat"))
NZDiffSub<-data_sum_NZ$ProLegStand[3] - data_sum_NZ$ProLegStand[2]
NZSESub<-data_sum_NZ$se[3]
NZ <-c("New Zealand", NZDiffSub, NZSESub, 25) #ADD quota threshold from Table 1

penaltyPro<-as.data.frame(rbind(Mexico, Peru, UK, Spain, USA, Port, Norway, Aus, France, Argentina, Brazil, NZ))

penaltySub$V2<-as.numeric(as.character(penaltySub$V2))
penaltySub$V3<-as.numeric(as.character(penaltySub$V3))
penaltySub$V4<-as.numeric(as.character(penaltySub$V4))

penaltyPro $V2<-as.numeric(as.character(penaltyPro $V2))
penaltyPro $V3<-as.numeric(as.character(penaltyPro $V3))
penaltyPro $V4<-as.numeric(as.character(penaltyPro $V4))

cor.test(penaltySub$V2, penaltySub$V4)
cor.test(penaltyPro $V2, penaltyPro $V4)

penaltySub
penaltyPro

##"Quota Penalty" - Difference between GBP and Q-GBP 

##Substantive Legitimacy 

plot3<-ggplot(penaltySub, aes(x= V4, y= V2, group= V1, label=V1)) + 
    geom_errorbar(aes(ymin= V2-V3, ymax= V2 + V3), width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, shape=21, fill="grey", aes(size=1), show.legend=FALSE) + # 21 is filled circle
    xlab("Quota Strength") +
   ylab("Treatment Effect Size (standardized w/in countries)") +
     ggtitle("Substantive Legitimacy 'Quota Penalty'") +
    theme_bw() +
    theme(legend.position=c(0.15,0.2))    +            # Position legend in bottom right
    scale_y_continuous(limits = c(-0.5, 0.25)) +
     scale_x_continuous(limits = c(0, 55)) +
  geom_hline(yintercept=0)+ 
  geom_label(size=2.5)


#Procedural Legitimacy 

plot4<-ggplot(penaltyPro, aes(x= V4, y= V2, group= V1, label=V1)) + 
    geom_errorbar(aes(ymin= V2-V3, ymax= V2 + V3), width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, shape=21, fill="grey", aes(size=1), show.legend=FALSE) + # 21 is filled circle
    xlab("Quota Strength") +
   ylab("Treatment Effect Size (standardized w/in countries)") +
     ggtitle("Procedural Legitimacy `Quota Penalty' ") +
    theme_bw() +
    theme(legend.position=c(0.15,0.2))    +            # Position legend in bottom right
    scale_y_continuous(limits = c(-0.5, 0.25))  +
       scale_x_continuous(limits = c(0, 55)) +
  geom_hline(yintercept=0)+ 
  geom_label(size=2.5)

require(gridExtra)

# Set the PDF output file
pdf("Figure3.pdf", width = 8, height = 5)  # Width and height are in inches

grid.arrange(plot3, plot4, ncol=2, nrow=1) 

# Close the graphics device
dev.off()



##FIGURE SI.1 
##Alternative measure

penaltySub $WomensRep <-recode(penaltySub $V1, "
'Argentina' = 45; 
'Australia'  = 38;   
'Brazil'    = 15;
'France'    = 37;
'Mexico'    = 48; 
'Norway'  =    44; 
'Peru'  = 40; 
'Portugal'     = 40; 
'Spain'  = 44 ; 
 'UK'       = 34; 
 'New Zealand' = 50;
 'USA' = 28")
 
penaltySub $AltMeasure<-penaltySub $V4 + penaltySub $WomensRep


penaltyPro $WomensRep <-recode(penaltyPro $V1, "
'Argentina' = 45; 
'Australia'  = 38;   
'Brazil'    = 15;
'France'    = 37;
'Mexico'    = 48; 
'Norway'  =    44; 
'Peru'  = 40; 
'Portugal'     = 40; 
'Spain'  = 44 ; 
 'UK'       = 34; 
 'New Zealand' = 50;
 'USA' = 28")
 
penaltyPro $AltMeasure<-penaltyPro $V4 + penaltyPro $WomensRep

##Substantive Legitimacy 

plot3<-ggplot(penaltySub, aes(x= AltMeasure, y= V2, group= V1, label=V1)) + 
    geom_errorbar(aes(ymin= V2-V3, ymax= V2 + V3), width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, shape=21, fill="grey", aes(size=1), show.legend=FALSE) + # 21 is filled circle
    xlab("Quota Strength") +
   ylab("Treatment Effect Size (standardized w/in countries)") +
     ggtitle("Substantive Legitimacy 'Quota Penalty'") +
    theme_bw() +
    theme(legend.position=c(0.15,0.2))    +            # Position legend in bottom right
    scale_y_continuous(limits = c(-0.5, 0.25)) +
     scale_x_continuous(limits = c(0, 100)) +
  geom_hline(yintercept=0)+ 
  geom_label(size=2.5)


#Procedural Legitimacy 

plot4<-ggplot(penaltyPro, aes(x= AltMeasure, y= V2, group= V1, label=V1)) + 
    geom_errorbar(aes(ymin= V2-V3, ymax= V2 + V3), width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, shape=21, fill="grey", aes(size=1), show.legend=FALSE) + # 21 is filled circle
    xlab("Quota Strength") +
   ylab("Treatment Effect Size (standardized w/in countries)") +
     ggtitle("Procedural Legitimacy `Quota Penalty' ") +
    theme_bw() +
    theme(legend.position=c(0.15,0.2))    +            # Position legend in bottom right
    scale_y_continuous(limits = c(-0.5, 0.25))  +
       scale_x_continuous(limits = c(0, 100)) +
  geom_hline(yintercept=0)+ 
  geom_label(size=2.5)


# Set the PDF output file
pdf("FigureSI1.pdf", width = 8, height = 5)  # Width and height are in inches

grid.arrange(plot3, plot4, ncol=2, nrow=1) 

# Close the graphics device
dev.off()


##SI Table 

## SI Table SI.7 

tab1<-penaltySub[,c(1, 2, 3, 4)]
tab2<-penaltyPro[,c(1, 2, 3, 4)]

tab<-cbind(tab1, tab2)

latex_table <- xtable(tab[,c(2, 3, 6, 7, 8)])
 

# Print LaTeX code to console
print(latex_table, type = "latex")

# Export to a .tex file
sink("tableSI.7.tex")
print(latex_table, type = "latex")
sink()

 
